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Abstract. In this short note, we present an easy to implement and fast algorithm for the 
computation of the steady solitary gravity wave solution of the free surface Euler equations in 
irrotational motion. First, the problem is reformulated in a fixed domain using the conformal 
mapping technique. Second, the problem is reduced to a single equation for the free surface. 
Third, this equation is solved using Petviashvili's iterations together with pseudo-spectral 
discretisation. This method has a super-linear complexity, since the most demanding opera- 
tions can be performed using a FFT algorithm. Moreover, when this algorithm is combined 
with the multi-precision arithmetics, the results can be obtained to any arbitrary accuracy. 



1. Introduction 



Solitary waves play a central role in nonlinear sciences. They appear in various fields 
ranging from plasmas physics to hydrodynamics and nonlinear optics [8, 16]. In some special 
r~| ■ cases, solitary wave solutions can be found analytically. For example, explicit expressions are 

known for integrable models such as KdV and NLS equations, but also for some non-integrable 
equations such as [13]. However, no closed form solutions are known for the practically very 
important case of the free surface Euler equations. In order to construct these solutions, 
one has to apply some approximation methods. Historically, asymptotic approximations have 
been proposed first [4] but, these series being divergent [7], the high-order approximations 
are only valid in the limit a ^ (a the wave amplitude). In order to avoid this limitation, 
numerical approaches have been proposed [5]. Several numerical methods have also been 
proposed [10]. One of the most widely used methods nowadays is Tanaka's algorithm [14]. 
In this note, we propose a simpler and faster method based on Petviashvili's iterations [12] 
leading to a very efficient numerical scheme for the computation of solitary gravity waves 
of the full Euler equations in the water of finite depth. Petviashvili's iterations have been 
successfully used for various simplified water wave equations but, to our knowledge, never 
for the fully nonlinear problem. The algorithm proposed here is efficient and it allows a very 
straightforward compact implementation under Matlab, for example. Such a script can be 
freely downloaded [1]. 

The note is organised as follows. In section 2, we introduce the hypothesis, the notations 
and the Babenko equation which is equivalent to the original problem. In section 3, we 
describe a very simple and fast algorithm to compute the solution, and we give some numerical 
results. Finally, the main conclusions are outlined in Section 4. 
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Figure 1. Sketch of physical and transformed domams. 



2. Mathematical model 



We consider a steady two-dimensional potential flow due to a solitary gravity wave prop- 
agating at the surface of a homogeneous fluid layer of constant depth d. Let be (x,y) a 
Cartesian coordinate system moving with the wave, x being the horizontal coordinate and 
y the upward vertical one. Since solitary waves are localised in space, the surface elevation 
tends to zero as x ^ ±cx3, and x = is the abscissa of the crest. The equations of the bottom, 
of the free surface and of the mean water level are given correspondingly hy y - -d, y - r/(x) 
and y - 0, respectively. 

2.1. Governing equation. Let be (j), ip, u and v the velocity potential, the stream function, 
the horizontal and vertical velocities, respectively, such that u = (px = ipy ^'^^ v - (py - -ipx- It 
is convenient to introduce the complex potential f = (p + iip (with i^ = -1) and the complex 
velocity w ^u-iv that are holomorphic functions z = x + iy (i.e., w - df/dz). The complex 
conjugate is denoted with a star (e.g., z* - x - \y), while subscripts 'b' denote the quantities 
written at the seabed — e.g., ZTa{x) - x - id, (ph{x) - (p{x,y--d) — and subscripts 's' denote 
the quantities written at the surface — e.g., Zs{x) - x + iT/(x), (ps{x) - (p{x,y = r]{x)). We 
emphasise that ips and V'b are constants because the surface and the bottom are streamlines. 

The far field velocity is such that (n, v) -> (-c, 0) as x ^ ±oo, so c is the wave phase velocity 
observed in the frame of reference where the fluid is at rest at infinity (c > if the wave travels 
to the increasing x-direction) . Note that c = (ipb - ips)/d due to the mass conservation. 

The dynamic condition can be expressed in term of the Bernoulli equation 



where p is the pressure divided by the density and g > is the acceleration due to gravity. At 
the free surface the pressure equals that of the atmosphere which is constant and set to zero 
without loss of generality, i.e., ps - 0- 

2.2. Conformal mapping. Let be the change of independent variable z C ^ (iips - f)/c, 
that conformally maps the fluid domain {-oo ^ x ^ oo;-d^y^ri} into the strip {-cxj ^ a ^ 
oo;-d ^ l3 ^ 0} where a = Re(C) and (3 = Im((^) (c.f. Fig. 1). The conformal mapping yields 
the Cauchy-Riemann relations Xa - yp and xp - -ya, while the complex velocity and the 
velocity components are 

(dzV^ -CXa - Cya 2 2 

Introducing new dependent variables X{a,j3) = x - a and Y{a,j3) = y - j3, the Cauchy- 
Riemann relations Xa - Yp and Xp - -Ya hold, while the bottom (/? = -d) and the free surface 



2p + 2gy + + V 




(2.1) 
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(/3 = 0) impermeabilities yield Y\y(a) e Y(a,~d) - and Ys{a) e Y(a,0) - rj. The functions 
X and Y can be expressed in term of Xy^ solely [3] 

X(a,/3) = cos[(/3 + d)a„]Xb(a), Y{a,/3) = sin[ (/3 + ] Xb(a), (2.2) 

so that the Cauchy-Riemann relations and the bottom impermeability are fulfilled identically. 
At the free surface (3-0, (2.2) yields Xs{a) = cos[dda] X\;,{a) , that can be inverted as 
Xh{a) - sec[dda] Xs{a), and hence the relation (2.2) yields 

Ys{a) = rj{a) = tan[ d9„ ] Xs(a), (2.3) 

which relates quantities written at the free surface only. The relation (2.3) can be trivially 
inverted giving, in particular, (daX)^ - ^{Y^} e daCot[dda]Ys, where is a self-adjoint 
pseudo-differential operator acting on a pure frequency as 

c^fjkax _ j kcoth{kd)e''''' (k^O) , . 

^^{i/d (k^o) 

This operator can be efficiently evaluated in the Fourier space using a FFT algorithm. 

2.3. Babenko's equation. A Lagrangian for water waves can be obtained considering the 
kinetic energy minus the potential one, leading, in the conformal plane, to the functional 
£ = f_°^ £ da with £ being the Lagrangian density [10] 

The governing equation is the Euler-Lagrange equation of this functional, i.e. 

° ^ If ^ '^{d^} ^ "'^^^^ ~ " ^^"^^^'^ " (^-^^ 

which is the Babenko equation for gravity solitary surface waves [2]. Applying the operator 
'ta^^ to (2.5) and splitting the linear and nonlinear parts, the equation can be presented as 

that is more convenient for the simple and efficient numerical resolution explained below. 

3. Numerical method and results 

First, we define the wave by choosing its Froude number F = c/\/gd. Second, we choose a 
KdV-like initial guess of the free surface elevation rjQ^a) 

r/o(a) = ^(F^- l)sech^(K;a/2), = tan{Kd) / Kd. (3.1) 

Third, equation (3.1-6) is used to compute k and thence to determine the length L of the 
computational domain (i.e. -L ^ a ^ L) such that rjQ^L) = to machine precision. The 
computational domain is then discretised with equally spaced nodes, with N large enough 
to avoid aliasing. Fourth, equation (2.6) is solved via Petviashvili's iterations [9, 11, 12] 

r^n.i - S^^-'{^{Vn}}, - fr^-^i^"!;!" , (3.2) 

J-c^Vnyy{vn}da 

where Sn is the co-called stabilisation factor which can be computed either in physical or 
Fourier spaces, thanks to the Parseval identity [15]. In the implementation, we systematically 
privilege the Fourier space since the operators ^ and can be very efficiently computed 
according to the definition (2.4) using any FFT algorithm. 
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The convergence of the iterative process (3.2) is checked by following the norm of the 
difference between two successive iterations along with the residual in the ioo norm 

where e„ are some prescribed tolerance parameters that are usually of the order of the floating 
point arithmetics precision. 

3.1. Numerical results. The present algorithm is illustrated on the following example of a 
fairly high amplitude solitary wave. It is validated and compared against Tanaka's algorithm 
[14]. Tanaka's solution is parametrized by the dimensionless parameter Qc - \uc\l\fgd e (0, 1), 
where Uc is the fluid velocity at the wave crest. We choose the special value qc - 0.4 for 
comparisons. Tanaka's algorithm was implemented in Matlab^ in its original version without 
any peculiar optimisations. The interval -30 ^ (/> ^ 30 is discretised using 20,001 nodes and 
the iterations are continued until the tolerance lO"^'^ is obtained between two successive values 
of the Froude number F, i.e. 

l-^n+i ~ \ < ^3 (- 10"^*^ in our computation). 

This computation required 47 iterations which lasted 385 s on our desktop computer. The 
resulting Froude number computed by Tanaka's algorithm is 

F = 1.274236297977903. 

Then, we take this Froude number and use it as the solitary wave parameter in the Babenko 
equation (2.6). We use a (periodised) domain of the same length as in Tanaka's method, and 
discretise it with = 16,384 nodes. The tolerance is chosen to be 10"^^, i.e. very close to 
machine accuracy. These parameters are kept in all the computations presented below. The 
iterative procedure is stopped when the loo norm of the difference between two successive 
iterations is less than the tolerance (the norm of the residual is checked a posteriori). Our 
algorithm required 372 iterations to converge and the computations lasted for about 0.5 s, i.e. 
almost three orders of magnitude faster than Tanaka's algorithm. This flagrant difference in 
CPU times can be explained by the fact that our method has the super-linear complexity 
0{N\ogN) while Tanaka's method has a quadratic complexity 0{N'^). 

Moreover, the present algorithm can be used to compute solitary waves to an arbitrary 
precision. To illustrate this, we implemented the scheme using a Multi-Precision Computing 
Toolbox (MPCT) for Matlab [6]. The code modifications needed to implement the arbi- 
trary accuracy are really minimalistic; this constitutes one of the major advantages of this 
MPCT. In these higher-precision experiments, all floating point operations are done with 50 
significant digits. The tolerance parameter is set to 10"^*^. To achieve such extreme accuracy, 
the Petviashvili scheme required 1,386 iterations. The convergence of the algorithm is shown 
on Figure 2-b. One can notice that the residual error is slightly higher than the difference 
between two successive iterations, likely due to the presence of quadratic nonlinearities into 
the equation. However, both measures of the numerical error decrease monotonically with 
iterations. The resulting solitary wave profile is depicted on Figure 2-a along with the cor- 
responding Tanaka solution. The two curves cannot be distinguished up to the graphical 
resolution, thus validating our method. The respective accuracies can be asserted comparing 
the dimensionless amplitudes ajd (Table 1). One can notice that despite the tolerance param- 
eter set to 10"^'^ in Tanaka's algorithm, its real accuracy is seriously degraded for this large 

"'^In the sequel, all the algorithms are compared in the same computing environment: Matlab®. 
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Figure 2. Solitary wave profile (a) and convergence of the numerical error during 
the iterations (b). 

Tanaka (16 digits) Present (16 digits) Present (50 digits) 

0.6819 448200954614 0.6819656673869021 0.68196566738691572070985777997577037303844037027343 
Table 1. Amplitude corresponding to F = 1.274236297977903 (wrong digits are underlined). 

amplitude wave (only five digits of the wave amplitude are predicted correctly). Conversely, 
the present method gives fourteen correct digits with a tolerance of 10^^^ in double precision 
calculation (about sixteen digits), and 49 correct digits using the 50-digit computation (the 
latter was determined using a 60-digit calculation with a tolerance of 10~^^ with eight times 
more nodes over an interval twice the size). 

The algorithm being validated, it can be used to produce some physically sound numerical 
results. For example, we compute to high-precision (here 10~^^, but it is not a limitation) 
the speed-amplitude relation for solitary waves. The result is shown on Figure 3. We also 
computed several values in a wide range of amplitudes using Tanaka's algorithm (represented 
with circles). Again, up to the graphical resolution we have a very good agreement between 
the two methods. However, the speed of Tanaka's algorithm does not allow to compute many 
values in a reasonable time. 

We note finally that, for simplicity, we defined the wave by its Froude number which varies 
in the interval 1 ^ F < 1.29421. Doing so, we can compute all waves such that ^ a < 0.796 
that is sufficient for all cases of practical interest. However, there are solutions of larger 
amplitudes, the limit a < 0.796 being obtained here because F is not a monotonic function 
of a (i.e., F{a) is multi-valued when a is large). Thus, in order to compute larger waves, 
another parameter (e.g., a/d or qc) should be used to define the wave. This requires simple 
modifications in the code [1]. We successfully did it and were able to compute all waves 
such that a < 0.82. For even higher waves, the present method is not efficient because these 
extreme waves being (almost) singular at the crest, the Fourier series does not converge fast 
enough and another spectral decomposition should be used instead. 
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Figure 3. Speed versus amplitude relation predicted by the present (solid line) and 
Tanaka's (circles) algorithms. 

4. Conclusions 

We have described a fast, accurate and easy to implement method for computing steady 
gravity sohtary surface waves solution of the full Euler equations. The method is based on first 
rewriting the problem into the equivalent but simpler Babenko's equation and, second, using 
Petviashvili's iterations to solve it numerically. The resulting scheme allows the computation 
to any arbitrary accuracy using a multi-precision floating-point arithmetics library [6]. All 
ingredients of the method are well known, but their combination seems to be new in this 
context and turns out to be very efficient. Our Matlab implementation is rather compact: 
The computational numerical core is about a dozen lines of code, the rest being pre- and post- 
processing. The script [1] can be freely downloaded and used (it can also be easily translated 
into any programming language). 
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